clear
clc
close all
%% 初始化
n = 1000;%迭代次数
fx = @(x) x(1)^2 + 2*x(2)^2 -2*x(1)*x(2);%函数
sx2 = @(x1) (2-x1^2)./x1;%边界
dfx = @(x) [2*x(1)-2*x(2);4*x(2)-2*x(1)];%梯度

H = [2,-2;-2,4];%Hessian矩阵

sample.x = [];
sample.dx = [];
sample.y = [];
x_iteration = repmat(sample,n,1);
%% 计算
x_iteration(1).x = [1;1];%起始点
x_iteration(1).y = fx(x_iteration(1).x);
x_iteration(1).dx = dfx(x_iteration(1).x);

s = -x_iteration(1).dx;%初始的sk
iteration = [1];
error = norm(x_iteration(1).dx);
Y = x_iteration(1).y;
figure
i = 1;
while(error(end)>=1e-30 && i<=n)
%for i = 2:1:n
    i = i+1;
    clf
    s0 = s;
    %梯度下降
    g = x_iteration(i-1).dx;
    alpha = -(g'*s)/(2*s'*H*s);
    x_iteration(i).x = x_iteration(i-1).x + alpha*s;

    %边界条件
    x_iteration(i).x(1) = min(x_iteration(i).x(1),sqrt(2)).*(x_iteration(i).x(1)>0);
    x_iteration(i).x(2) = ((x_iteration(i).x(2) < sx2(x_iteration(i).x(1)))...
        .*(x_iteration(i).x(2) -sx2(x_iteration(i).x(1)))+sx2(x_iteration(i).x(1))).*(x_iteration(i).x(2)>0);

    x_iteration(i).y = fx(x_iteration(i).x);

    %求梯度
    x_iteration(i).dx = dfx(x_iteration(i).x);
    beta = (norm(x_iteration(i).dx)^2)/(norm(x_iteration(i-1).dx)^2);
    s = -x_iteration(i).dx+beta*s0;

    %绘图
    iteration = [iteration,i];
    error = [error,norm(x_iteration(i).dx)];
    Y = [Y,x_iteration(i).y];
    subplot(2,1,1)
    semilogy(iteration,error);
    yticks(10.^(-30:3:0));
    subplot(2,1,2)
    semilogy(iteration,Y);
    yticks(10.^(-60:5:0));
    pause(0.001);
end